#delimit;
clear;
version 12;
set more off, perm;
  quietly log;
  local logon = r(status);
  if "`logon'" == "on" {; log close; };
log using SS_SMPO_epa.log, text replace;

/*******************************/
/*  Mark David Nieman  */
/*  1/21/2016          */
/*  Application: EPA data */
/*  Stata 12.1        */
/*******************************/

#delimit;
clear;
clear matrix;
clear mata;
program drop _all;
set seed 1999;


/* strategic logit with partial observability */
	#delimit;
	program define slpo_lf, rclass;
		args lnf A_L A_Rr B_Rr;
		tempvar pB_r pA_R ;
			quietly gen double `pB_r' = (exp(`B_Rr'))/(1+exp(`B_Rr'));
			quietly gen double `pA_R' = (exp(`pB_r'*`A_Rr'))/(exp(`A_L')+exp(`pB_r'*`A_Rr'));			
		quietly replace `lnf' = ln((1-`pA_R') + (`pA_R')*(1-`pB_r')) if $ML_y1==0;
		quietly replace `lnf' = ln((`pA_R')*(`pB_r')) if $ML_y1==1;
	end;
	
/*load data*/
use CAA_Dataset1.dta, clear;
/* set=1 if pp or hospital*/

sort afsid year;
gen id=_n if afsid!=afsid[_n-1];
replace id=id[_n-1] if afsid==afsid[_n-1];
drop if id==.;
drop if afsid==afsid[_n-1] & year==year[_n-1];
isid id year;
tsset id year;

/* logit analysis */
logit nc_year major  public PercentBlack PercentHisp PercentPoverty MedHsInc PercentUnemploy L.nc_year;

matrix ll_l = e(ll);	
gen sample = 1 if e(sample)==1;
sum nc_year if sample==1;
tab nc_year if sample==1;
sutex nc_year major  public PercentBlack PercentHisp PercentPoverty MedHsInc PercentUnemploy if sample==1, minmax nobs;

egen N = count(sample) if sample==1;
egen Obs = max(N);
drop N;
			/* Model fit */
			/* i's likelihood */
			gen double lnfj_l = . ;
			replace lnfj_l = ln(invlogit(-(_b[major]*major +  _b[public]*public + _b[PercentBlack]*PercentBlack +
										   _b[PercentHisp]*PercentHisp + _b[PercentPoverty]*PercentPoverty + _b[MedHsInc]*MedHsInc + 
										   _b[PercentUnemploy]*PercentUnemploy + _b[L.nc_year]*L.nc_year + _b[_cons]))) if nc_year==0 & sample==1;
			
			replace lnfj_l = ln(invlogit( (_b[major]*major +  _b[public]*public + _b[PercentBlack]*PercentBlack +
										   _b[PercentHisp]*PercentHisp + _b[PercentPoverty]*PercentPoverty + _b[MedHsInc]*MedHsInc + 
										   _b[PercentUnemploy]*PercentUnemploy + _b[L.nc_year]*L.nc_year + _b[_cons]))) if nc_year==1 & sample==1;

			gen double pr_vio_l = invlogit(_b[major]*major  + _b[public]*public + _b[PercentBlack]*PercentBlack +
										   _b[PercentHisp]*PercentHisp + _b[PercentPoverty]*PercentPoverty + _b[MedHsInc]*MedHsInc + 
										   _b[PercentUnemploy]*PercentUnemploy + _b[L.nc_year]*L.nc_year + _b[_cons]) if sample==1;
			
			gen pred_vio_l = 0 if sample==1;
				replace pred_vio_l =1 if pr_vio_l > .5 & sample==1;
			gen corr_l = 0 if sample==1;
				replace corr_l = 1 if (pred_vio_l==1 & nc_year==1) | (pred_vio_l==0 & nc_year==0);
			gen corr_vio_l = .;
				replace corr_vio_l = 0 if (pred_vio_l==0 & nc_year==1);
				replace corr_vio_l = 1 if (pred_vio_l==1 & nc_year==1);
			gen incorr_vio_l =.;
				replace incorr_vio_l = 0 if (pred_vio_l==0 & nc_year==0);
				replace incorr_vio_l = 1 if (pred_vio_l==1 & nc_year==0);		
			sum corr_l corr_vio_l incorr_vio_l;							   
											   
#delimit;
/* get initial conditions */
logit nc_year public PercentBlack PercentHisp PercentPoverty MedHsInc PercentUnemploy ;
	matrix B2 = e(b);
	predict xb;
	gen L = exp(xb)/(1+exp(xb));

	gen Lmajor =-major;
	gen Lpp = -pp;	
	gen R_pub = L*public;
	gen R_pun = L*L.pun_t;
	gen R_lag = L*L.nc_year;
	gen R_ins = L*L.insp_t;
	gen R_black = L*PercentBlack;
	gen R_hisp = L*PercentHisp;
	gen R_pov = L*PercentPoverty;
	gen R_inc = L*MedHsInc;
	gen R_unemp = L*PercentUnemploy;
	gen R_c = L;
#delimit;
logit nc_year Lmajor R_black R_hisp R_pov R_inc R_unemp R_lag R_c , nocons;
	matrix B1 = e(b);

/*matrices of yr and state fes*/
matrix yr = (0,0,0,0,0,0,0,0,0,0,0,0);
matrix st = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);

#delimit;	
/* Strategic logit with partial observability */
ml model lf slpo_lf (A_L: nc_year = major , nocons) 
					(A_Rr:  = PercentBlack PercentHisp PercentPoverty MedHsInc PercentUnemploy L.nc_year) 
					(B_Rr:  = public PercentBlack PercentHisp PercentPoverty MedHsInc PercentUnemploy ) , tech(nr) ;
ml init B1 B2, copy ;
ml maximize, diff;	
matrix ll_slpo = e(ll);
gen sample_m = 1 if e(sample)==1;
			
			/* Model fit */
			/* i's likelihood */
			gen double xb_A_L = ([A_L]_b[major]*major );
			gen double xb_A_Rr = ([A_Rr]_b[PercentBlack]*PercentBlack + [A_Rr]_b[PercentHisp]*PercentHisp + [A_Rr]_b[PercentPoverty]*PercentPoverty +
										[A_Rr]_b[MedHsInc]*MedHsInc + [A_Rr]_b[PercentUnemploy]*PercentUnemploy + [A_Rr]_b[L.nc_year]*L.nc_year + [A_Rr]_b[_cons]);
			gen double xb_B_Rr = ([B_Rr]_b[public]*public + [B_Rr]_b[PercentBlack]*PercentBlack + [B_Rr]_b[PercentHisp]*PercentHisp + [B_Rr]_b[PercentPoverty]*PercentPoverty +
										[B_Rr]_b[MedHsInc]*MedHsInc + [B_Rr]_b[PercentUnemploy]*PercentUnemploy + [B_Rr]_b[_cons]);
			gen double A_L = exp(xb_A_L)/(exp(xb_A_L)+exp(xb_A_Rr));
			gen double A_Rr = exp(xb_A_Rr)/(exp(xb_A_L)+exp(xb_A_Rr));
			gen double B_Rr = invlogit(xb_B_Rr);
			
			
			gen double lnfj_slpo = .;
			replace lnfj_slpo = ln(A_L + A_Rr*(1-B_Rr)) if nc_year == 0 & sample_m==1;
			replace lnfj_slpo = ln(A_Rr*B_Rr) if nc_year == 1 & sample_m==1;
			
			gen pr_san_slpo = A_Rr*B_Rr;
			gen pr_det_slpo = A_L;
			gen pr_acq_slpo = A_Rr*(1-B_Rr);
			gen pred_vio_slpo = 0 if sample_m==1;
				replace pred_vio_slpo =1 if pr_san_slpo>max(pr_det_slpo,pr_acq_slpo) & sample_m==1;
			gen corr_slpo = 0 if sample_m==1;
				replace corr_slpo = 1 if (pred_vio_slpo==1 & nc_year==1) | (pred_vio_slpo==0 & nc_year==0);
			gen corr_vio_slpo = .;
				replace corr_vio_slpo = 0 if (pred_vio_slpo==0 & nc_year==1);
				replace corr_vio_slpo = 1 if (pred_vio_slpo==1 & nc_year==1);
			gen incorr_vio_slpo =.;
				replace incorr_vio_slpo = 0 if (pred_vio_slpo==0 & nc_year==0);
				replace incorr_vio_slpo = 1 if (pred_vio_slpo==1 & nc_year==0);		
			sum corr_slpo corr_vio_slpo incorr_vio_slpo;	
			
			gen pred_Y3_slpo = A_Rr*(1-B_Rr);
			gen Acq_gt_33 = 0 if sample_m==1;
				replace Acq_gt_33 = 1 if pr_acq_slpo>0.30 & sample_m==1;
			tab Acq_gt_33;
			

			/* Voung Test*/
			svmat ll_l;
			svmat ll_slpo;
			gen num_V = ll_slpo1 - ll_l1 + ((15/2)*ln(Obs)-(9/2)*ln(Obs));
			egen denom1_V_1 = mean((lnfj_slpo/lnfj_l)^2);
			egen denom1_V_2 = mean(lnfj_slpo/lnfj_l);
			gen denom1_V_3 = denom1_V_2^2;
			gen denom2_V = sqrt(denom1_V_1 - denom1_V_3);
			gen Voung = num_V/denom2_V;
			sum num_V denom2_V Voung;

			
			/* Clarke test */
			signtest lnfj_slpo=lnfj_l;		
			
			/* Summary fit stats */
			sum corr_slpo corr_vio_slpo incorr_vio_slpo;
			sum corr_l corr_vio_l incorr_vio_l;
			sum pred_Y3_slpo;
			tab Acq_gt_33;
	
	log off;
